/** \file funcoes.cpp
 *  \brief Implementação das funções definidas no arquivo funcoes.hpp
 *  \author Geraldo Francisco de S. Rebouças | gfs.reboucas@gmail.com
 *  \since 01/12/2012
 *  \version 1.0
 *  \date 02/12/2012 */

#include "funcoes.hpp"

void solver(double T0, double TM, vetor X0, double W, double PV, vetor (*F)(double, vetor, double))
{/** \param T0 Recebe o valor inicial da variável independente, normalmente o tempo
  *  \param TM Recebe o valor final da variável independente
  *  \param X0 Recebe o vetor de estados iniciais
  *  \param W Recebe a frequência da excitação senoidal e com ele calcula o passo de integração
  *  \param PV Recebe o parâmetro que desejamos variar durante as simulações, muito útil quando traçamos o diagrama de bifurcação do sistema
  *  \param F Recebe a função que representa a EDO a ser resolvida    */

    double hh, h, t;
    hh = 2.0*M_PI/W;
    h = hh/1000.0;

//    scout << h << "\n";

    int c = 0; // contador para a seção de Poincaré

    ofstream rst("RST.txt"); // Declaração e associação dos arquivos de texto aos seus respectivos objetos.
    ofstream pfs("PFS.txt");
    ofstream pcr("PCR.txt"); /// para que PCR.txt renda uma boa imagem, deverá ter 1Mb de tamanho, ou 3*1e4 pontos. Isso equivale a k*1e3*(3*1e4) iterações
    ofstream opi("OPI-2.txt");

    t = T0; // Condições iniciais
    vetor X(X0);

    while(t <= TM)
    {
        if(t >= 0.5*TM) /// k = 0.5 (vide comentario do arquivo PCR.txt)
        {
            pfs << scientific << t << "\t" << X << "\n"; // Dica: deixar comentado p/ TM grande ou comportamento caótico
//            if(c % 1000 == 0)
//                pcr << scientific << t << "\t" << X << "\n";
        }
        rst << scientific << t << "\t" << X << "\n"; // Dica: Deixar comentado p/ TM grande

//        if((t >= 2.577050) && (t <= 2.577053))
//            opi << scientific << t << "\t" << X << "\n";

        RungeKutta4(h, t, X, PV, F);
        c++;
    }
    rst.close(); // Fechamento dos arquivos de texto
    pfs.close();
    pcr.close();

//    for(signed int i = 1; i < 5; i++)
//    OPI("PCR.txt", i);
}

void solverVE(double T0, double TM, vetor X0, double W, double PV, vetor (*F1)(double, vetor, double), // Eq. do modelo
                                                                       matriz (*F2)(double, vetor, matriz, double)) // Eq. Variacional
{/** \param T0 Recebe o valor inicial da variável independente, normalmente o tempo
  *  \param TM Recebe o valor final da variável independente
  *  \param X0 Recebe o vetor de estados iniciais
  *  \param W Recebe a frequência da excitação senoidal e com ele calcula o passo de integração
  *  \param PV Recebe o parâmetro que desejamos variar durante as simulações, muito útil quando traçamos o diagrama de bifurcação do sistema
  *  \param F1 Recebe a função que representa a EDO associada a equação variacional
  *  \param F2 Recebe a função que representa a Equação variacional    */

    double h, t;
    h = 2.0*M_PI/W;
    h = h/1000.0;

    ofstream rst("RST.txt"); // Declaração e associação dos arquivos de texto aos seus respectivos objetos.

    t = T0; // Condições iniciais
    vetor X(X0);
    matriz M;
    M.setIdent(X0.retDim());

    while(t <= TM)
    {
        rst << scientific << t << "\t" << X << M << "\n";
        RungeKutta4VE(h, t, X, M, PV, F1, F2);
    }

    rst.close();
}

void solverLyap(double T0, double TM, vetor X0, double W, double PV, vetor (*F1)(double, vetor, double), // Eq. do modelo
                                                                       matriz (*F2)(double, vetor, matriz, double)) // Eq. Variacional
{/** \param T0 Recebe o valor inicial da variável independente, normalmente o tempo
  *  \param TM Recebe o valor final da variável independente
  *  \param X0 Recebe o vetor de estados iniciais
  *  \param W Recebe a frequência da excitação senoidal e com ele calcula o passo de integração
  *  \param PV Recebe o parâmetro que desejamos variar durante as simulações, muito útil quando traçamos o diagrama de bifurcação do sistema
  *  \param F1 Recebe a função que representa a EDO associada a equação variacional
  *  \param F2 Recebe a função que representa a Equação variacional    */

    double h, t;
    h = 2.0*M_PI/W;
    h = h/1000.0;

    int q = X0.retDim();

    ofstream rst("RST.txt"); // Declaração e associação dos arquivos de texto aos seus respectivos objetos.
    ofstream lpn("LPN.txt");

    t = T0; // Condições iniciais
    vetor X(X0);
    matriz M;
    M.setIdent(q);

    matriz U(q, q);
    U.setIdent(q);
    matriz V(q, q);
    vetor lmb(q); // Lambda
    vetor s(q); // Somatório
    matriz dx(q, q);
    int k = 0;

    vetor tmpv(q);
    double nm;

    while(t <= TM)
    {
        k++;
        dx = M*U;
        for(int i = 0; i < q; i++)
        {
            V.setCol(dx.col(i), i); // V(:,i) = dx(:,i)
            for(int j = 0; j < i; j++)
            {
                nm = V.col(i)^U.col(j);
                tmpv = V.col(i) - nm*U.col(j);
                V.setCol(tmpv, i);
            }
            tmpv = V.col(i);
            nm = tmpv.norma(2);
            tmpv = V.col(i)*(1.0/nm);
            U.setCol(tmpv, i);
            s(i) = s(i) + log(nm);
            lmb(i) = s(i)/(h*k);
        }
        rst << scientific << t << "\t" << X << M << "\n";
        lpn << scientific << t << "\t" << lmb << "\n";

        RungeKutta4VE(h, t, X, M, PV, F1, F2);
    }

    rst.close();
    lpn.close();
    lmb.show();
}

void solverDBF(double T0, double TM, vetor X0, double W, double PI, double PF, int S, vetor (*F)(double, vetor, double))
{/** \param T0 Recebe o valor inicial da variável independente, normalmente o tempo
  *  \param TM Recebe o valor final da variável independente
  *  \param X0 Recebe o vetor de estados iniciais
  *  \param W Recebe a frequência da excitação senoidal e com ele calcula o passo de integração
  *  \param PI Recebe o valor inicial parâmetro que desejamos variar durante as simulações
  *  \param PF Recebe o valor final parâmetro que desejamos variar durante as simulações
  *  \param S Recebe a quantidade de iterações a serem simuladas variando PV
  *  \param F Recebe a função que representa a EDO a ser resolvida    */

    double h, t;
    h = 2*M_PI/W;
    h = h/1000.0;
    int c; // Contador para o diagrama de Bifurcação
    double hp = (PF - PI)/S; // passo do parâmetro variável
    ofstream dbf("DBF.txt");

    while(PI <= PF)
    {
        t = T0;
        vetor X(X0);
        c = 0;
        while(t <= TM)
        {
            if((t >= 0.8*TM) && (c % 1000 == 0))
                dbf << scientific << PI << "\t" << X << "\n";

            RungeKutta4(h, t, X, PI, F);
            c++;
        }
        PI += hp;
        dbf << endl;
    }
    dbf.close();
}

void RungeKutta4(double H, double &T, vetor &X, double PV, vetor (*F)(double, vetor, double))
{/** \param H Recebe o valor do passo de integração
  *  \param T Recebe o valor da variável independente
  *  \param X Recebe o vetor de estados
  *  \param PV Recebe o parâmetro que desejamos variar durante as simulações, muito útil quando traçamos o diagrama de bifurcação do sistema
  *  \param F Recebe a função que representa a EDO a ser resolvida    */
    int o = X.retDim();
    vetor k1(o), k2(o), k3(o), k4(o);

    k1 = H*F(T, X, PV);
    k2 = H*F(T + 0.5*H, X + 0.5*k1, PV);
    k3 = H*F(T + 0.5*H, X + 0.5*k2, PV);
    k4 = H*F(T + H, X + k3, PV);

    T += H;
    X += (k1 + 2.0*(k2 + k3) + k4)*(1.0/6.0);
}

void RungeKutta4VE(double H, double &T, vetor &X, matriz &M, double PV, vetor (*F1)(double, vetor, double), // Eq. do modelo
                                                                           matriz (*F2)(double, vetor, matriz, double)) // Eq. Variacional
{
    int dd = X.retDim();
    int ll = M.retl();
    int cc = M.retc();

    vetor k1(dd), k2(dd), k3(dd), k4(dd);
    matriz K1(ll, cc), K2(ll, cc), K3(ll, cc), K4(ll, cc);


    k1 = H*F1(T, X, PV);
    K1 = H*F2(T, X, M, PV);

    k2 = H*F1(T + 0.5*H, X + 0.5*k1, PV);
    K2 = H*F2(T + 0.5*H, X + 0.5*k1, M + 0.5*K1, PV);

    k3 = H*F1(T + 0.5*H, X + 0.5*k2, PV);
    K3 = H*F2(T + 0.5*H, X + 0.5*k2, M + 0.5*K2, PV);

    k4 = H*F1(T + H, X + k3, PV);
    K4 = H*F2(T + H, X + k3, M + K1, PV);

    T += H;
    X += (k1 + 2.0*(k2 + k3) + k4)*(1.0/6.0);
    M += (K1 + 2.0*(K2 + K3) + K4)*(1.0/6.0);
}


void OPI(const char *pc, int k)
{/** \param pc Recebe o arquivo de texto contendo os pontos da seção de Poincaré
  *  \param k Recebe o valor do período desejado    */
    int tam = qtLin(pc);

    ifstream pcr(pc);
    ofstream dst("DST.txt");

    vector <double> t(tam);
    vector <double> x1(tam);
    vector <double> x2(tam);

    vector <double> d(tam - k);
    vector <double> t1(tam - k);
    vector <double> t2(tam - k);

    double a, b;

    for(int i = 0; i < tam; i++) // leitura dos pontos da seção de Poincaré
        pcr >> t[i] >> x1[i] >> x2[i];

    pcr.close();

    for(int j = 0; j < tam - k; j++)
    {
        a = pow(x1[j] - x1[j + k], 2);
        b = pow(x2[j] - x2[j + k], 2);

        dst << scientific << sqrt(a + b) << "\t" << t[j] << "\t" << t[j + k] << "\n";
    }
    dst.close();

    ifstream dist2("DST.txt");

    for(int l = 0; l < tam - k; l++) // leitura dos pontos do arquivo  DST.txt
        dist2 >> d[l] >> t1[l] >> t2[l];

    kMenor(d, t1, t2, k);
}

void kMenor(vector <double> A, vector <double> B, vector <double> C, int k)
{/// http://blog.teamleadnet.com/2012/07/quick-select-algorithm-find-kth-element.html - Adam Horvath

    if(A.size() < k)
    {
        cerr << "Erro! Dimensoes incompativeis! \n";
        abort();
    }

    int a = 1;
    int b = A.size() - 1;
    int c;
    int d;
    double e;
    double fa, fb, fc;

    while(a < b) /// se a == b alcançamos o k-ésimo elemento
    {
        c = a;        d = b;        e = A[0.5*(c + d)];

        while(c < d) /// para se c e d coincidirem
        {
            if(A[c] >= e) /// coloca o maior valor no fim
            {
                fa = A[d];
                A[d] = A[c]; /// parte que eu preciso ajeitar p o algoritmo das OPIS
                A[c] = fa;

                fb = B[d];
                B[d] = B[c]; /// parte que eu preciso ajeitar p o algoritmo das OPIS
                B[c] = fb;

                fc = C[d];
                C[d] = C[c]; /// parte que eu preciso ajeitar p o algoritmo das OPIS
                C[c] = fc;

                d--;
            }else /// o valor é maior do que o pivô, não importa
                c++;
        }

        if(A[c] > e) /// se passamos direto (c++) precisamos retornar
            c--;

        if(k <= c) /// o ponteiro c está no final dos primeiros k elementos
            b = c;
        else
            a = c + 1;
    }
    scout << k <<"a menor distancia: "<< A[k] << " u. c. \n";
    scout << "Instante inicial: " << B[k] << " s \t Instante final: " << C[k] << " s \n";
}

int qtLin(const char *tmp)
{/** \param tmp Recebe o nome do arquivo de texto   */
    int C, T = 0;
    FILE* A = fopen(tmp, "r");

    while((C = fgetc(A)) != EOF)
    {
        if(C == '\n')
            T++;
    }

    return T;
}

double Rand(double a, double b)     {    srand (time(NULL));    return (b - a)*(rand() / double(RAND_MAX)) + a;     }
/** \fn Rand
 *  \brief Gera um número pseudo-aleatório no intervalo \f$ (a,b) \f$
 *  \param a Menor elemento do conjunto dos possíveis números que podem ser gerados
 *  \param b Maior elemento do conjunto dos possíveis números que podem ser gerados */

void beep(int a)
{/** \param a Recebe a quantidade de tempo que o beep durará    */
    clock_t cfinal;
    cfinal = clock() + a*CLOCKS_PER_SEC;
    while(clock() < cfinal)
        cout << "\a";
}
